home *** CD-ROM | disk | FTP | other *** search
/ Languguage OS 2 / Languguage OS II Version 10-94 (Knowledge Media)(1994).ISO / gnu / gmp-132.lha / gmp-1.3.2 / mpz_gcd.c < prev    next >
C/C++ Source or Header  |  1993-05-19  |  4KB  |  170 lines

  1. /* mpz_gcd -- Calculate the greatest common divisior of two integers.
  2.  
  3. Copyright (C) 1991 Free Software Foundation, Inc.
  4.  
  5. This file is part of the GNU MP Library.
  6.  
  7. The GNU MP Library is free software; you can redistribute it and/or modify
  8. it under the terms of the GNU General Public License as published by
  9. the Free Software Foundation; either version 2, or (at your option)
  10. any later version.
  11.  
  12. The GNU MP Library is distributed in the hope that it will be useful,
  13. but WITHOUT ANY WARRANTY; without even the implied warranty of
  14. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15. GNU General Public License for more details.
  16.  
  17. You should have received a copy of the GNU General Public License
  18. along with the GNU MP Library; see the file COPYING.  If not, write to
  19. the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.  */
  20.  
  21. #include "gmp.h"
  22. #include "gmp-impl.h"
  23. #include "longlong.h"
  24.  
  25. #ifndef BERKELEY_MP
  26. void
  27. #ifdef __STDC__
  28. mpz_gcd (MP_INT *w, const MP_INT *u, const MP_INT *v)
  29. #else
  30. mpz_gcd (w, u, v)
  31.      MP_INT *w;
  32.      const MP_INT *u;
  33.      const MP_INT *v;
  34. #endif
  35. #else /* BERKELEY_MP */
  36. void
  37. #ifdef __STDC__
  38. gcd (const MP_INT *u, const MP_INT *v, MP_INT *w)
  39. #else
  40. gcd (u, v, w)
  41.      const MP_INT *u;
  42.      const MP_INT *v;
  43.      MP_INT *w;
  44. #endif
  45. #endif /* BERKELEY_MP */
  46. {
  47.   mp_size usize, vsize, wsize;
  48.   mp_ptr up_in, vp_in;
  49.   mp_ptr up, vp;
  50.   mp_ptr wp;
  51.   mp_size i;
  52.   mp_limb d;
  53.   int bcnt;
  54.   mp_size w_bcnt;
  55.   mp_limb cy_digit;
  56.  
  57.   usize = ABS (u->size);
  58.   vsize = ABS (v->size);
  59.  
  60.   /* GCD(0,v) == v.  */
  61.   if (usize == 0)
  62.     {
  63.       if (w->alloc < vsize)
  64.     _mpz_realloc (w, vsize);
  65.  
  66.       w->size = vsize;
  67.       MPN_COPY (w->d, v->d, vsize);
  68.       return;
  69.     }
  70.  
  71.   /* GCD(0,u) == u.  */
  72.   if (vsize == 0)
  73.     {
  74.       if (w->alloc < usize)
  75.     _mpz_realloc (w, usize);
  76.  
  77.       w->size = usize;
  78.       MPN_COPY (w->d, u->d, usize);
  79.       return;
  80.     }
  81.  
  82.   /* Make U odd by shifting it down as many bit positions as there
  83.      are zero bits.  Put the result in temporary space.  */
  84.   up = (mp_ptr) alloca (usize * BYTES_PER_MP_LIMB);
  85.   up_in = u->d;
  86.   for (i = 0; (d = up_in[i]) == 0; i++)
  87.     ;
  88.   count_leading_zeros (bcnt, d & -d);
  89.   bcnt = BITS_PER_MP_LIMB - 1 - bcnt;
  90.   usize = mpn_rshift (up, up_in + i, usize - i, bcnt);
  91.  
  92.   bcnt += i * BITS_PER_MP_LIMB;
  93.   w_bcnt = bcnt;
  94.  
  95.   /* Make V odd by shifting it down as many bit positions as there
  96.      are zero bits.  Put the result in temporary space.  */
  97.   vp = (mp_ptr) alloca (vsize * BYTES_PER_MP_LIMB);
  98.   vp_in = v->d;
  99.   for (i = 0; (d = vp_in[i]) == 0; i++)
  100.     ;
  101.   count_leading_zeros (bcnt, d & -d);
  102.   bcnt = BITS_PER_MP_LIMB - 1 - bcnt;
  103.   vsize = mpn_rshift (vp, vp_in + i, vsize - i, bcnt);
  104.  
  105.   /* W_BCNT is set to the minimum of the number of zero bits in U and V.
  106.      Thus it represents the number of common 2 factors.  */
  107.   bcnt += i * BITS_PER_MP_LIMB;
  108.   if (bcnt < w_bcnt)
  109.     w_bcnt = bcnt;
  110.  
  111.   for (;;)
  112.     {
  113.       int cmp;
  114.  
  115.       cmp = usize - vsize != 0 ? usize - vsize : mpn_cmp (up, vp, usize);
  116.  
  117.       /* If U and V have become equal, we have found the GCD.  */
  118.       if (cmp == 0)
  119.     break;
  120.  
  121.       if (cmp > 0)
  122.     {
  123.       /* Replace U by (U - V) >> cnt, with cnt being the least value
  124.          making U odd again.  */
  125.  
  126.       usize += mpn_sub (up, up, usize, vp, vsize);
  127.       for (i = 0; (d = up[i]) == 0; i++)
  128.         ;
  129.       count_leading_zeros (bcnt, d & -d);
  130.       bcnt = BITS_PER_MP_LIMB - 1 - bcnt;
  131.       usize = mpn_rshift (up, up + i, usize - i, bcnt);
  132.     }
  133.       else
  134.     {
  135.       /* Replace V by (V - U) >> cnt, with cnt being the least value
  136.          making V odd again.  */
  137.  
  138.       vsize += mpn_sub (vp, vp, vsize, up, usize);
  139.       for (i = 0; (d = vp[i]) == 0; i++)
  140.         ;
  141.       count_leading_zeros (bcnt, d & -d);
  142.       bcnt = BITS_PER_MP_LIMB - 1 - bcnt;
  143.       vsize = mpn_rshift (vp, vp + i, vsize - i, bcnt);
  144.     }
  145.     }
  146.  
  147.   /* GCD(U_IN, V_IN) now is U * 2**W_BCNT.  */
  148.  
  149.   wsize = usize + w_bcnt / BITS_PER_MP_LIMB + 1;
  150.   if (w->alloc < wsize)
  151.     _mpz_realloc (w, wsize);
  152.  
  153.   wp = w->d;
  154.  
  155.   MPN_ZERO (wp, w_bcnt / BITS_PER_MP_LIMB);
  156.  
  157.   cy_digit = mpn_lshift (wp + w_bcnt / BITS_PER_MP_LIMB, up, usize,
  158.               w_bcnt % BITS_PER_MP_LIMB);
  159.   wsize = usize + w_bcnt / BITS_PER_MP_LIMB;
  160.   if (cy_digit != 0)
  161.     {
  162.       wp[wsize] = cy_digit;
  163.       wsize++;
  164.     }
  165.  
  166.   w->size = wsize;
  167.  
  168.   alloca (0);
  169. }
  170.